#############################################################
########################### SET-UP ##########################
#############################################################

# load packages
Packages <- c('ggplot2', 'tidyverse', 'scales', 'DiagrammeR', 'ggdag', 
              'knitr', 'pander', 'foreign', 'dplyr', 'ggeffects', 'car', 
              'gridExtra', 'grid', 'texreg', 'viridis', 'stargazer', 
              'nnet', 'margins', 'emmeans', 'effects', 'mediation', 
              'broom', 'ggridges', 'patchwork', 'fuzzySim', 'ggrepel', 
              'effectsize') 

pacman::p_load(Packages, character.only = TRUE)

# load data
load(file = "Replication_dataset.rda")

#############################################################
#################### FIGURES 1 and 2 ########################
#############################################################

# Recode
ag_viz <- opinion_model %>% 
  mutate(agegroup = ifelse(age >= 20 & age <= 29, "20-29",
                                         ifelse(age >= 30 & age <= 39, "30-39",
                                                       ifelse(age >= 40 & age <= 49, "40-49",
                                                                     ifelse(age >= 50  & age <= 59, "50-59",
                                                                                   ifelse(age >= 60 & age <= 69, "60-69",
                                                                                                 ifelse(age >= 70 & age <= 79, "70-79",
                                                                                                               ifelse(age >= 80 & age <= 89, "80-89", NA)))))))) %>% 
  dplyr::select(cohort, year, age, 
                prefsep_educ,  prefsep_hc,  prefsep_childcare,  prefsep_elderly,  
                prefsep_energy,  prefsep_env, prefsep_defence,  prefsep_justice,  
                prefsep_arts,  prefsep_regions,  prefsep_trans, prefsep_farmers, 
                prefsep_welfare,  prefsep_jobprogram,  prefsep_sspoor) %>% 
  mutate(id = row_number(),
         id = as.numeric(id)) %>% 
  gather(policy, response, prefsep_educ:prefsep_sspoor) %>% 
  mutate(want_more = ifelse(response == 1, 1, 
                            ifelse(response == -1 | response == 0, 0, NA)),
         want_less = ifelse(response == 1 | response == 0, 0, 
                            ifelse(response == -1, 1, NA)),
         want_sq = ifelse(response == -1 | response == 1, 0, 
                          ifelse(response == 0, 1, NA))) %>% 
  group_by(age, year) %>% 
  mutate(mean_more = mean(want_more, na.rm = T),
         mean_less = mean(want_less, na.rm = T),
         mean_sq = mean(want_sq, na.rm = T)) %>% 
  ungroup() %>% 
  distinct(age, year, mean_more, mean_less, mean_sq) %>% 
  drop_na(age) %>% 
  gather(type, value, mean_more:mean_sq) %>% 
  mutate(type = ifelse(type == "mean_more", "Wants more spending", 
                       ifelse(type == "mean_less",  "Wants less spending", "Wants the same amount of spending")))

# Figure 1
figure1 <- ggplot(aes(x = age, y = value), data = subset(ag_viz,type == "Wants more spending")) +
  geom_point(alpha = .1) +
  geom_smooth(color = "black")+
  theme_minimal() +
  scale_y_continuous(limits = c(.3,.6),
                     labels = scales::percent_format(accuracy = 1))+
  labs(x = "Age",
       y = "Share of respondents") +
  facet_wrap(~year)

figure1 

# Figure 2
figure2 <- ggplot(aes(x = age, y = value), data = subset(ag_viz,type == "Wants the same amount of spending")) +
  geom_point(alpha = .1) +
  geom_smooth(color = "black")+
  theme_minimal() +
  scale_y_continuous(limits = c(.3,.6),
                     labels = scales::percent_format(accuracy = 1))+
  labs(x = "Age",
       y = "Share of respondents") +
  facet_wrap(~year)

figure2

#############################################################
######################### FIGURE 3 ##########################
#############################################################

# Recode - Support for more spending
ag_more <- opinion_model %>% 
  dplyr::select(cohort, age, year, norm_income, gender, generation,
                prefsep_educ,  prefsep_hc,  prefsep_childcare,  prefsep_elderly,  prefsep_energy,  prefsep_env,  
                prefsep_defence,  prefsep_justice,  prefsep_arts,  prefsep_regions,  prefsep_trans,  
                prefsep_farmers,  prefsep_welfare,  prefsep_jobprogram,  prefsep_sspoor, religiosity, party, 
                norm_educ, employment, marital) %>% 
  mutate(id = row_number(),
         id = as.numeric(id)) %>% 
  gather(policy, response, prefsep_educ:prefsep_sspoor) %>% 
  mutate(response = ifelse(response == 0, 0,
                           ifelse(response == -1, 0, 
                                  ifelse(is.na(response) == T, NA, 1)))) %>% 
  group_by(id) %>% 
  mutate(prefmore_spending = mean(response, na.rm = T)) %>% 
  ungroup() %>% 
  distinct(id, age, generation, year, prefmore_spending, norm_income, gender, cohort, religiosity, party, norm_educ, employment,marital) %>% 
  mutate(year_f = as.factor(year),
         party = ifelse(party == "bloc", "bq", party))

# Model without control for year and generation
mod1 <- lm(prefmore_spending ~ age + I(age^2) + norm_income + gender + religiosity + party + norm_educ + employment + marital, data = ag_more)

# Model with control for year and generation
mod2 <- lm(prefmore_spending ~ age + I(age^2) + norm_income + gender + year + generation + religiosity + party + norm_educ + employment + marital, data = ag_more)

# Predicted values
p_age <- ggeffect(mod1, terms = c("age [25:75]")) %>% mutate(model = "Age only")
p_age2 <- ggeffect(mod2, terms = c("age [25:75]")) %>% mutate(model = "Controlling for generations and years") %>%  full_join(p_age) 

# Figure 3 - More spending
p<-ggplot(aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, linetype = model), data = p_age2) +
  geom_line(size = 1) +
  geom_ribbon(alpha = .1) +
  scale_y_continuous(limits = c(.35,.55),
                     labels = scales::percent_format(accuracy = 1))+
  scale_linetype_manual(name = "", values = c("dotted", "solid"))+
  labs(subtitle = "Spend more",
       x = "Age",
       y = "") +
  theme_minimal()+
  theme(legend.position = "")

p

# Recode - Support for the same amount of spending
ag_sq <- opinion_model %>% 
  dplyr::select(cohort, age, year, norm_income, gender, generation,
                prefsep_educ,  prefsep_hc,  prefsep_childcare,  prefsep_elderly,  prefsep_energy,  prefsep_env,  
                prefsep_defence,  prefsep_justice,  prefsep_arts,  prefsep_regions,  prefsep_trans,  
                prefsep_farmers,  prefsep_welfare,  prefsep_jobprogram,  prefsep_sspoor, religiosity, party, norm_educ, employment,marital) %>% 
  mutate(id = row_number(),
         id = as.numeric(id)) %>% 
  gather(policy, response, prefsep_educ:prefsep_sspoor) %>% 
  mutate(response = ifelse(response == 1, 0,
                           ifelse(response == -1, 0, 
                                  ifelse(is.na(response) == T, NA, 1)))) %>% 
  group_by(id) %>% 
  mutate(prefsq_spending = mean(response, na.rm = T)) %>% 
  ungroup() %>% 
  distinct(id, age, generation, year, prefsq_spending, norm_income, gender, cohort, religiosity, party, norm_educ, employment,marital) %>% 
  mutate(year_f = as.factor(year),
         party = ifelse(party == "bloc", "bq", party))

# Model without control for year and generation
mod3 <- lm(prefsq_spending ~ age + I(age^2) + norm_income + gender + religiosity + party + norm_educ + employment + marital, data = ag_sq)

# Model with control for year and generation
mod4 <- lm(prefsq_spending ~ age + I(age^2) + norm_income + gender + year + generation + religiosity + party + norm_educ + employment + marital, data = ag_sq)

# Predicted values
p_age <- ggeffect(mod3, terms = c("age [25:75]")) %>% mutate(model = "Age only")
p_age2 <- ggeffect(mod4, terms = c("age [25:75]")) %>% mutate(model = "Controlling for\ngenerations\nand years") %>%  full_join(p_age)

# Figure 3 - Same amount of spending
p2<-ggplot(aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, linetype = model), data = p_age2) +
  geom_line(size = 1) +
  geom_ribbon(alpha = .1) +
  scale_y_continuous(limits = c(.35,.55),
                     labels = scales::percent_format(accuracy = 1))+
  scale_linetype_manual(name = "", values = c("dotted", "solid"))+
  labs(subtitle = "Spend the same amount",
       x = "Age",
       y = "") +
  theme_minimal()+
  theme(legend.position = "right")

p2

# Patch plots together and save
figure3 <- p + p2 
figure3

#############################################################
####################### FIGURES 4-5-6 #######################
#############################################################

# Recode
models <- opinion_model %>% 
  dplyr::select(cohort, age, year, norm_income, gender, generation,religiosity, party, norm_educ, employment,marital,
                prefsep_educ,  prefsep_hc,  prefsep_childcare,  prefsep_elderly,  prefsep_energy,  prefsep_env,  
                prefsep_defence,  prefsep_justice,  prefsep_arts,  prefsep_regions,  prefsep_trans,  
                prefsep_farmers,  prefsep_welfare,  prefsep_jobprogram,  prefsep_sspoor) %>% 
  mutate(id = row_number(),
         id = as.numeric(id)) %>% 
  gather(policy, response, prefsep_educ:prefsep_sspoor) %>% 
  mutate(prefless = ifelse(response == 1, 0,
                           ifelse(response == 0, 0, 
                                  ifelse(is.na(response) == T, NA, 1))),
         prefmore = ifelse(response == -1, 0,
                           ifelse(response == 0, 0, 
                                  ifelse(is.na(response) == T, NA, 1))),
         prefsq = ifelse(response == 1, 0,
                         ifelse(response == -1, 0, 
                                ifelse(is.na(response) == T, NA, 1))),
         party = ifelse(party == "bloc", "bq", party)) %>% 
  dplyr::select(-id, -response, -cohort) %>% 
  arrange(policy)

# Estimate models - Support for more spending or the same level of spending - 15 policies separately
names <- unique(models$policy)

estmore <- models %>%
  group_by(policy) %>%
  group_map(~ lm(prefmore ~ age + I(age^2) + generation + year + norm_income + gender + party + religiosity + employment + norm_educ + marital, data = .)) 
names(estmore) <- names
estsq <- models %>%
  group_by(policy) %>%
  group_map(~ lm(prefsq ~ age + I(age^2) + generation + year + norm_income + gender + party + religiosity + employment + norm_educ + marital, data = .))
names(estsq) <- names

# Predicted values from all models
listmore <- list()
listsq <- list()
for (i in 1:15) {
  listmore[[i]] <- ggeffect(estmore[[i]], terms = c("age [25:75]"))
  listsq[[i]] <- ggeffect(estsq[[i]], terms = c("age [25:75]"))
}
names(listmore) <- names
names(listsq) <- names

prefmore <- bind_rows(listmore, .id="df") %>% 
  as.data.frame() %>% 
  mutate(df = str_remove_all(df, "prefsep_"),
         group = "spending more")
prefall <- bind_rows(listsq, .id="df") %>% 
  as.data.frame()%>% 
  mutate(df = str_remove_all(df, "prefsep_"),
         group = "spending the same amount") %>% 
  full_join(prefmore)  %>% 
  mutate(df = ifelse(df == "arts", "Arts",
                     ifelse(df == "childcare", "Childcare",
                            ifelse(df == "defence", "Defence", 
                                   ifelse(df == "educ", "Education",
                                          ifelse(df == "elderly", "Serv. for the elderly", 
                                                 ifelse(df == "energy", "Energy", 
                                                        ifelse(df == "env", "Environment",  
                                                               ifelse(df == "farmers", "Farmers",  
                                                                      ifelse(df == "hc", "Healthcare",
                                                                             ifelse(df == "jobprogram", "Job-creation prog.",
                                                                                    ifelse(df == "justice", "Justice", 
                                                                                           ifelse(df == "regions", "Regions", 
                                                                                                  ifelse(df == "sspoor", "Social serv. for the poor", 
                                                                                                         ifelse(df == "trans", "Transportation",
                                                                                                                ifelse(df == "welfare", "Welfare", NA))))))))))))))),
         group = factor(group, levels = c("spending the same amount", "spending more")))

# Figure 4
figure4 <- ggplot(aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, linetype = group), data = subset(prefall, str_detect(df, "Justice|Transportation|Defence|elderly|Healthcare"))) + 
  geom_line(size = 1) +
  geom_ribbon(alpha = .2) +
  scale_x_continuous(limits = c(25,75), breaks = c(25,50,75)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_linetype(name = "Government should be...") +
  theme_minimal() +
  labs(x = "Age", 
       y = "") +
  facet_wrap(~df, ncol = 3) +
  theme(legend.position = "bottom")

figure4

# Figure 5
figure5a <- ggplot(aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, linetype = group), data = subset(prefall, str_detect(df, "Education|Childcare|Job"))) + 
  geom_line(size = 1) +
  geom_ribbon(alpha = .2) +
  scale_x_continuous(limits = c(25,75), breaks = c(25,50,75)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_linetype(name = "Government should be...") +
  theme_minimal() +
  labs(x = "Age", 
       y = "") +
  facet_wrap(~df, ncol = 3) +
  theme(legend.position = "bottom")

figure5b <- ggplot(aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, linetype = group), data = subset(prefall, str_detect(df, "Environment"))) + 
  geom_line(size = 1) +
  geom_ribbon(alpha = .2) +
  scale_x_continuous(limits = c(25,75), breaks = c(25,50,75)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_linetype(name = "Government shoulde be...") +
  theme_minimal() +
  labs(x = "Age", 
       y = "") +
  facet_wrap(~df, ncol = 1) +
  theme(legend.position = "")

figure5 <- figure5a + figure5b + 
  plot_layout(widths = c(3, 1))
figure5

# Figure 6
figure6 <- ggplot(aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, linetype = group), data = subset(prefall, str_detect(df, "Regions|poor|Farmers|Energy|Welfare|Arts"))) + 
  geom_line(size = 1) +
  geom_ribbon(alpha = .2) +
  scale_x_continuous(limits = c(25,75), breaks = c(25,50,75)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_linetype(name = "Government should be...") +
  theme_minimal() +
  labs(x = "Age", 
       y = "") +
  facet_wrap(~df, ncol = 3) +
  theme(legend.position = "bottom")

figure6 

#############################################################
########################## TABLE 2 ##########################
#############################################################

listmore2 <- list()
listless2 <- list()
listsq2 <- list()
for (i in 1:15) {
  listmore2[[i]] <- ggeffect(estmore[[i]], terms = c("age [20:70]", "generation"))
  listsq2[[i]] <- ggeffect(estsq[[i]], terms = c("age [20:70]", "generation"))
}
names(listmore2) <- names
names(listsq2) <- names

prefmore2<- bind_rows(listmore2, .id="df") %>% 
  as.data.frame() %>% 
  mutate(df = str_remove_all(df, "prefsep_"),
         policy = "spending more")

prefall2 <- bind_rows(listsq2, .id="df") %>% 
  as.data.frame()%>% 
  mutate(df = str_remove_all(df, "prefsep_"),
         policy = "spending the same amount") %>% 
  full_join(prefmore2)  %>% 
  mutate(df = ifelse(df == "arts", "Arts",
                     ifelse(df == "childcare", "Childcare",
                            ifelse(df == "defence", "Defence", 
                                   ifelse(df == "educ", "Education",
                                          ifelse(df == "elderly", "Serv. for the elderly", 
                                                 ifelse(df == "energy", "Energy", 
                                                        ifelse(df == "env", "Environment",  
                                                               ifelse(df == "farmers", "Farmers",  
                                                                      ifelse(df == "hc", "Healthcare",
                                                                             ifelse(df == "jobprogram", "Job-creation prog.",
                                                                                    ifelse(df == "justice", "Justice", 
                                                                                           ifelse(df == "regions", "Regions", 
                                                                                                  ifelse(df == "sspoor", "Social serv. for the poor", 
                                                                                                         ifelse(df == "trans", "Transportation",
                                                                                                                ifelse(df == "welfare", "Welfare", NA))))))))))))))),
         policy = factor(policy, levels = c("spending the same amount", "spending more")))  

scenarios <- prefall2 %>% 
  filter(str_detect(df, "Education|Healthcare|elderly|Transportation|Environment")) %>% 
  filter((group == "3third" & x >= 32 & x <= 91) | (group == "4fourth" & x >= 18 & x <= 63) ) %>% 
  filter(x == 20 | x == 30 | x == 40 | x == 50 | x == 60 | x == 70) %>% 
  filter(policy == "spending more") %>% 
  dplyr::select(df, group, x, predicted) %>% 
  mutate(predicted = predicted*100) %>% 
  spread(group, predicted) 

# 1928-1955 birth cohort
diffs_3third <- scenarios %>% 
  dplyr::select(-`4fourth`) %>% 
  filter(x == 30 | x == 40 | x == 50 | x == 60 | x == 70) %>% 
  replace(is.na(.), 0) %>% 
  group_by(df) %>%
  mutate(`3third` = `3third` - lag(`3third`, default = first(`3third`), order_by = x)) %>% 
  as.data.frame() 

diffs_3third

# 1928-1955 birth cohort
diffs_4fourth <- scenarios %>% 
  dplyr::select(-`3third`) %>% 
  filter(x == 20 | x == 30 | x == 40 | x == 50 | x == 60) %>% 
  replace(is.na(.), 0) %>% 
  group_by(df) %>%
  mutate(`4fourth` = ifelse(x!=20,`4fourth` - lag(`4fourth`, default = first(`4fourth`), order_by = x), `4fourth`))  %>% 
  as.data.frame()

diffs_4fourth
